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ABSTRACT 

We explore a variety of statistics of clusters selected with cosmic shear measurement 
by utilizing both analytic models and large numerical simulations. We first develop 
a halo model to predict the abundance and the clustering of weak lensing selected 
clusters. Observational effects such as galaxy shape noise are included in our model. 
We then generate realistic mock weak lensing catalogs to test the accuracy of our 
analytic model. To this end, we perform full-sky ray-tracing simulations that allow 
us to have multiple realizations of a large continuous area. We model the masked 
regions on the sky using the actual positions of bright stars, and generate 200 mock 
weak lensing catalogs with sky coverage of ^1000 squared degrees. We show that 
our theoretical model agrees well with the ensemble average of statistics and their 
covariances calculated directly from the mock catalogues. With a typical selection 
threshold, ignoring shape noise correction causes overestimation of the clustering of 
weak lensing selected clusters with a level of about 10%, and shape noise correction 
boosts the cluster abundance by a factor of a few. We calculate the cross-covariances 
using the halo model with accounting for the effective reduction of the survey area 
due to masks. The covariance of the cosmic shear auto power spectrum is affected by 
the mode-coupling effect that originates from sky masking. Our model and the results 
can be readily used for cosmological analysis with ongoing and future weak lensing 
surveys. 
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1 INTRODUCTION 


The accelerating expansion of the universe is now established by an array of astronomical observations such as Type la 
supernovae (e.g., [Betoule et al. 2014[), measurement of baryon acoustic oscillations in galaxy surveys (e.g. 


2011 


2013 


Blake et al. 12011 


Planck Collaboration et al. 


Anderson et al.|2014|, anisotropies of the cosmic microwave background (CMB) (e.g., 


Beutler et al. 


Hinshaw et al. 


20141, large-scale galaxy distribution (e.g., Sanchez et al. | |2012| [Beutler et al.||2014[ ), and 


weak gravitational lensing (e.g., Kilbinger et al.|2013 |. In order to realize the cosmic acceleration within the theory of general 
relativity, an exotic form of energy needs to be postulated to dominate in the present-day universe. There is another possibility 
to explain the cosmic acceleration without dark energy, e.g., modified gravity theory. Modihed gravity models do not assume 
an unknown energy but change essentially the basic equation of gravitational action. Observationally, measurement of the 
growth of matter density fluctuations will help us to distinguish the models including the Einstein gravity with dark energy, 
because the modification of gravity induce characteristic clustering patterns in matter density distribution. 
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2 M. Shirasaki et al. 


Gravitational leasing is a powerful probe of the matter distribution in the universe. Small image distortions of distant 
galaxies are caused by intervening mass distribution. Small distortion caused by the large-scale structure of the universe is 
called cosmic shear. It contains, in principle, rich information on the matter distribution at small and large scales and the 
evolution over time. Image distortion induced by gravitational lensing is, however, very small in general. Therefore, we need 
statistical analyses of the cosmic shear signal by sampling a large number of distant galaxies in order to extract cosmological 
information from gravitational lensing. The conventional statistic of cosmic shear is two-point correlation function or its 
Fourier-counterpart, power spectrum. If the cosmic shear field obeys a Gaussian distribution, the two-point statistics suffice 
to describe all the information of cosmic shear. However, this is not the case in reality because cosmic shear has non-Gaussian 
information caused by non-linear gravitational growth ( Sato et al.||2009 K In order to extract the full information content, it 
is desirable to use other statistical quantities that probes nonlinear structure of length scale of ~10 Mpc or less. Glusters of 
galaxies are one the most reliable objects for this purpose. 

The number count of clusters is expected to be highly sensitive to growth of matter density perturbations ( Lilje|[T9^ l, 
whereas the spatial correlation of the position of clusters and cosmic shear provides the information on the matter density 
profile as well as clustering of clusters (e.g., Oguri et al. 2012 Okabe et al. 2013 Covone et al.|[20l4 l. Fortunately, cosmic 
shear itself provides an efficient way of locating galaxies of clusters. Cluster finding methods with cosmic shear are based on 
reconstruction of matter density distribution over an area of sky (Hamana, Takada fc Yoshida|200^ Hennawi fc Spergel|2005 


Maturi et al. 2005 Marian et al.|[2012K A reconstructed mass density map can be used to identify high density regions as 


“peaks” that are mostly caused by massive collapsed objects such as clusters of galaxies ( [Miyazaki et al.|2007[ jSchirmer et al.j 
2007 Shan et al.|2012 1. The unique advantage of weak lensing among various techniques is that it does not rely on uncertain 


physical state of the baryonic component in clusters. There have been a number of studies that investigate cosmological 
information in number counts of weak lensing selected clusters (Maturi et al.|20T0 Kratochvil, Haiman fc May|2010| Dietrich 


fc Hartla^|2010 Yang et al.||2011 Hilbert et al.|[20l2 l. Recently, Marian et al. (20131 combined other statistics beyond the 
abundance of weak lensing selected clusters. The authors conclude that correlation analysis of weak lensing selected clusters 
allow one to derive tight constraints on cosmological parameters. 

In the present paper, we study in detail the properties of a class of statistics of weak lensing selected clusters. Our study 
is aimed at being applied to real observations such as Subaru Hyper-Suprime-Cam Survey. It is well-known that the intrinsic 
ellipticities of source galaxies induce noise to lensing shear maps. The so-called shape noise causes typically false detection of 
clusters with cosmic shear measurement. We develop theoretical framework to model the correspondence of underlying dark 
matter halos and weak lensing selected clusters in presence of shape noise. Sky masking causes another important observational 
effect on statistical analyses of weak lensing selected clusters ( |Liu et ^|2014b| |. Since reconstructed mass density is usually 
defined by local cosmic shear signals, boundaries of masked regions would make the reconstruction inaccurate. In order to 
realize the realistic situation in galaxy imaging surveys, we perform gravitational lensing simulations on curved full-sky. We 
then utilize these simulations to create two hundreds of mock weak lensing catalogs with the proposed sky coverage in ongoing 
Hyper Suprime-Gam (HSC) survejQ For a large set of mock HSC surveys, we generate reconstructed mass map and identify 
the local maxima on each map as an indicator of weak lensing selected clusters. These simulations enable us to study the 
statistical property of weak lensing selected clusters in presence of shape noise and masked region. We are also able to examine 
our theoretical model through the large set of realistic mock weak lensing observations. 

The rest of the paper is organized as follows. In Section we describe the methodology to search clusters with cosmic 
shear measurement. There, we present the statistical property of weak lensing selected clusters and the theoretical model 
of statistics of interest. In Section]^ we use a large set of A^-body simulations to perform full-sky lensing simulations and 
create mock weak lensing maps incorporated with the information of HSC surveys. In Section we provide the result of 
our measurement of statistical quantities over a set of full-sky and masked sky simulations. We also compare the simulation 
results and our theoretical models in detail. Conclusions and discussions are summarized in Section 0 


2 WEAK LENSING 

We summarize the basics of weak gravitational lensing effect in this section. We also describe the finder algorithm of galaxy 
clusters with weak lensing measurement. 


1 
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Probing cosmology with weak tensing selected clusters I 3 


2.1 Basics 


When considering the observed position of a source object as G and the true position as (3, one can characterize the distortion 
of image of a source object by the following 2D matrix: 


— = ( 1 - - 7i -72 \ 

86^ ~ y —72 1 — K + 71 y ’ 


( 1 ) 


where k is convergence and 7 is shear. 

One can relate each component of Aij to the second derivative of the gravitational potential as follows (Bartelmann & 
|Schneider|200l)[Munshi et al.|2008| ); 
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where x is the comoving distance and r(x) represents the comoving angular diameter distance. Gravitational potential <E> can 
be related to matter density perturbation 5 according to Poisson equation. Therefore, convergence can be expressed as the 
weighted integral of 5 along the line of sight; 




^X9{X,X)-- 


( 5 ) 


2.2 Cluster finding 


Weak lensing provides a physical method to reconstruct the projected matter density field. The conventional technique for 
reconstruction is based on the smoothed map of cosmic shear. Let us first define the smoothed convergence field as 


k:{ 0 ) = J <P9' K{e - e')uie'), 


( 6 ) 


where U is the filter function to be specified below. Although we can calculate the same quantity by smoothing the shear field 
7 (e.g., Shirasaki fc Yoshida|2014 1, we use Eq. § for simplicity in the following. 

Various functional form of U are proposed in literature (e.g., Hamana, Takada fc Yoshida|2004 Hennawi fc Spergel|2005 


Maturi et al.|20lo Marian et al.|2012). We consider the Gaussian filteij^ 




( 7 ) 


With this filter, we can easily model the statistical properties of the contaminant of a smoothed K. map, called shape noise. 
The noise in a Af map would follow the Gaussian distribution when one can use a sufficient large number of source galaxies 
and when source galaxies are oriented randomly. The Gaussian properties of the noise makes it easy to model the lensing 
peak statistics, as will be shown in the following. 

We denote the shape noise contribution to a smoothed lensing map by A/”. For a given smoothing scale 9g, correlation 
function of the shape noise after Gaussian smoothing is given by ( [van Waerbeke|2000[ ) 

| 6 »- 6 »'P' 


mew{0')) = 




exp 


29% 


( 8 ) 


where a-y is the rms of the intrinsic ellipticity of sources and rigai represents the number density of source galaxies. One can 
derive the power spectrum of noise convergence field JV by Fourier transforming of Eq. 0; 


PAfit) = 


2ngai 


exp 




Using Eq. @, we define the moment of Af as 


''noise,I 






1/2 


( 9 ) 


( 10 ) 


^ In practice, the filter function should be compensated in order to remove an undetermined constant convergence < |Schneider|l996| l. In 
Appendix]^ we examine compensated Gaussian filters when searching for weak-lensing clusters. There, we show that our model can be 
suitably modified for the case of compensated Gaussian filters. 
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In a smoothed leasing map, peaks with high signal-to-noise ratio v = /C/crnoise,o are likely associated with cluster of 
galaxies (e.g., Hamana, Takada fc Yoshida||2004 1. We first locate high peaks on a /C map and then relate each peak with an 
isolated massive halo along the line of sight. We assume the following universal density profile of dark matter halos (Navarro, 
Frenk fc White|1997 1: 

= V , .2' ^ (11) 


{r/rs){l + r/rsf' 

where Vs and ps represent the scale radius and the scale density, respectively. The parameters and pa can be essen¬ 
tially convolved into one parameter, the concentration Cvii{M,z), by the use of two halo mass relations; namely, M = 
47 rrvijAvir(i;)pcrit(2)/3, where rvir is the virial radius corresponding to the overdensity criterion Avir(2) as shown in, e.g., 


Navarro, Frenk & White (19971, and M — j dV ph{rs,pa) with the integral performed out to rvir. In this paper, we adopt the 


functional form of the concentration parameter in Duffy et al. (20081 

- 0.081 


Cvir(M, 2 ) = 5.72 


M 


( 1 O 14 / 1 - 1 M 0 


( 1 +^) 


- 0.71 


The corresponding convergence can be calculated as in Hamana, Takada Yoshida (20041, 

(P\- ^psrafiR/ra) 

•^crit 

where R represents the perpendicular proper distance from the center of halo and f{x) is 
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In Eq. ( [T^ , Ecrit is defined by the following relation 
"" 4)i^ A As ’ 

where Dg, D\, and Dig are the angular diameter distance to the source, to the lens, and between the source and the lens, 
respectively. 

In order to predict the peak height in tC map, we need to take the following effects into account: (i) the offset between 
the position of a peak and the center of the corresponding halo and (ii) the modulation of peak height due to the shape noise. 


Fan, Shan & Liu (20101 have studied these two effects using numerical simulations and analytic approach. Let us hrst work 


on the simple assumption that the peak position is set to be the halo center. The peak height in absence of shape noise is 
given by 


^peak,lt J" d eu{6-,6G)Kh{0). 


(16) 


The actual peak height on a noisy IC map is not given by Eq. (161, but it obeys a probability distribution (Fan, Shan & Liu 


|2010p . The probability distribution function for a given /Cpeak.ii is calculated by 

T-, 1 /iy~ \ \ ^peak.N I A.'peak.obs A.'peak,?! j 

Prob(A^peak,obs|A-peak,/i) = 


(17) 


/ f*peak,N (^peak.obs I ^Peak,/i )dAip|,g^jj 

where /Cpeak.oba is the measured peak height and Upeak.N is defined as the expected number density of peaks with the measured 
peak height of /Cpeak.oba when the halo contribution /Cpeak.h is known in advance. For derivation of ripgak.N, we decompose the 
observed peak height /Cpeak.oba into three components: 


A3peak,obs — M -|- A^LSS + ^peak.li, 


(18) 


where N is the noise convergence field caused by shape noise, /Clss and /Cpeak.h represent the convergence field due to 
large-scale structure and foreground halos, respectively. Note that /Cpeak.h is a known quantity to derive ripeak.N- We aim at 
determining the relationship between /Cpeak.obs and a given A^peak,h(z, Af). 


Following Fan, Shan & Liu (20101, we assume that the noise field N is given by a Gaussian distribution with the power 
spectrum of Eq. If /Clss is a Gaussian random field, at the position of peaks, the total noise field (i.e. JV + IC-lss) obeys the 
probability distribution function of Gaussian peaks. Also, we can calculate the contribution from the (known) corresponding 
halo once the difference between the peak position and the halo center is specified. We assume that the peak position is at 
the center of the corresponding halo, although the equality does not hold in general. We have checked that the assumption is 
indeed reasonable for peaks with high signal-to-noise ratio in the case of Oq ~ 2arcmin, cr-y = 0.4, and Ugai ~ 10arcmin“^. 









































Probing cosmology with weak tensing selected clusters I 5 



Figure 1. The effective selection function of galaxy clusters. Black line shows the expected convergence signal caused by an isolated 
dark matter halo. Red line shows the modulated peak height on a noisy smoothed convergence map (see the text for the definition). We 
plot each line in units of the variance of shape noise cTnoise o = 0.017. 


Therefore, we set npeak,N to be the number density of peaks for Gaussian field A/* + /Clss- Using the relation of A/* + /Clss = 
A^peak,obs — A^peak,h, we Can obtain the number density npeak,N as (see, |Fan, Shan &; Liu|p010[ l for details) 
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(19) 


where Ci is the ith moment of A/” + /Clss, = 2((Ti/(J2)^, 'fN = rrJ/((ToO' 2 ) and /Cpeak.h denotes the second derivative of 
ICpea,k,h with respect to 9i at the halo centre. Here, F{xN\ICpeak,h) is calculated as follows: 


F lx N\Klpeak,h') — exp 


K. 
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In Eqs. (191 and (201, xm and ejv are given by 
Aivi + Ajv 2 Aivi + Ajv2 


xn = 


ejv 


( 20 ) 


( 21 ) 


(72 2a2Xi^ 

where Xni represents the ith diagonal component of the second derivative tensor of A/” + /Clss- 

We adopt a-y — 0.4, rigai = 10arcmin“^ and the source redshift is set to be ^source = 1. These are the typical values for 
ground-based galaxy imaging surveys (e.g., Heymans et al.|2012 |. Also, we employ a Gaussian smoothing with the full width 
at half maximum of 5 arcmin. This corresponds to 9a = 5/\/81n2 = 2.12 arcmin and thus to (Tnoise.o — 0.017. Although the 
smoothing scale adopted here is slightly larger than that in previous works (see, e.g., Hamana, Takada fc Yoshida||2004 1 by 
a factor of about two, the noise level on the smoothed convergence map is similar. Gaussian smoothing with ~ 2 arcmin 
with the actual data set is already examined in |Shan et al.|p012[ ). 

Let us examine the effect of shape noise on weak lensing peaks. For this purpose, we dehne the mean modulation of peak 
height in a noisy fC map as follows: 


Afpeak.obs)^, M) = / d/C/C Prob(/C|/Cpeak,h(2,M)). 


( 22 ) 


Figurej^shows the comparison with ^peak,obs(2, M) and ICpeak,hlz, M) for a given dark matter halo with mass of M at redshift 
z. In this figure, red line shows the contour of ^peak,obs( 2 , M) in units of (Jnoiae.o, whereas black line indicates the contour of 
tCpeak,hlz, M). In a noisy K. map, the shape noise modulates the height of peaks and the number of peaks slightly increases. 
We have tested the validity of our model against numerical simulations. The result is shown in Appendix [X] 

As shown in Figure[^ the Gaussian smoothing of ~ 2 arcmin are effective to search for clusters with mass of ~ Mq 
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at 2 ~ 0.1 — 0.2. The selection of mass and redshift is basically determined by the typical angular size of dark matter halos of 
interest (e.g., Hamana, Takada fc Yoshida||2004 1 . Naively, it is expected that lower redshift clusters are detected with larger 
smoothing scales. In order to verify this expectation, we have studied the statistical properties of leasing peaks when adopting 
a Gaussian smoothing with the full width at half maximum of 15 arcmin (corresponding to Oq ~ 6.4 arcmin). In this case, we 
do not find one-to-one correspondence between selected peaks and halos, and thus our analytic model does not work. This is 
likely caused by the so-called projection effect; the effective redshift of leasing in the case of ^source = 1 is 0.1 — 0.5 whereas we 
attempt to search for clusters at lower redshift 2 < 0.1. We argue that our model is valid when the smoothing scale is set to 
be 1 — 2 arcmin, corresponding to the typical angular size of dark matter halos with mass of ~ Mq at 2 ~ 0.1 — 0.5. 


2.3 Statistics 

We consider a set of statistics derived from weak leasing measurement. In order to extract cosmological information from 
the number and the distribution of massive dark matter halos, we utilize peaks on a smoothed leasing map as described in 
Section 12.21 


Convergence power spectrum 


First, we consider the power spectrum of convergence. Under the flat sky approximation, the Fourier transform of convergence 
field is defined by 

k(0) = 


(2 




(23) 


The power spectrum of convergence field is defined by 

{k{li)k{t2)) = (24) 

where is the Dirac delta function. By using Limber approximatioij^ ( Limber||l954 Kaiser|[l992 1 and Eq. ([^, one can 

calculate the convergence power spectrum as follows: 


p.4i) = 


I 


WJyf 
dx r.i Ps 




(25) 


r(x)2 

where Ps{k) is the three dimensional matter power spectrum, Xs is comoving distance to source galaxies and IFk(x) is the 
leasing weight function defined as 
2 


W4x) = 


Ho 

c 


II mO 


r{Xs - x)r{x) 

r{Xs) 


(i + ^(x)). 


(26) 


The non-linear gravitational growth of density fluctuations significantly affects the amplitude of convergence power 
spectrum at the angular scales less than 1 degree ( Jain, Seljak fc White|[2000 Hilbert et al.|[2009 Sato et 5!||2009 K Typical 
weak leasing surveys aim at measuring the angular scales larger than a few arcmin corresponding to a few Mpc. This leads 
weak leasing can be one of the most powerful probes for constraints on dark matter distribution at Mpc scale. Therefore, 
accurate theoretical prediction of non-linear matter power spectrum is essential for deriving cosmological constraints from 
weak lensing power spectrum. In order to predict the non-linear evolution of Ps {k) for standard ACDM universe, a numerical 


approach based on Wbody simulations has given steady results over the past few decades (Peacock & Dodds 1996 


et al. 


et al. 


2003 


( 20121 . 


Heitmann et al.||2010 Takahashi et al.||2012l. We adopt the most recent model of non-linear Ps{k) of 


Smith 


Takahashi 


Convergence peak count 

The number count of massive clusters is sensitive to various cosmological parameters such as the equation of state of dark 
energy (e.g., Allen, Evrard fc Mantz|[2011 1 . In this section, we use peak counts as a cosmological probe. We locate the local 
maxima in a smoothed lensing map and associate each identified peak with a massive dark matter halo along the same line 
of sight. 

In practice, one can select a lensing peak by its peak height. We define the signal-to-noise ratio of a peak hj v = 


® The validity of Limber approximation have been discussed in e.g., |jeong, Komatsu fc Jain| |2009[ |. The typical accuracy of Limber 
approximaion is of a level of 1% for t > 10. 
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A^peak,obs/o'noise,o- For a given threshold r'thre, one can predict the surface number density of peaks with v > r-thre as follows 
(e.g., Hamana, Takada fc Yoshida||2004 1: 




POO 

J O’ 


d./Cpeak,obs Prob(/Cpeak,obs I AIpeak,/i (2^5 1 


(27) 


^thre ^noise,0 


where dn/dM represents the mass function of dark matter halo and the volume element is expressed as (PV/AzAQ, = /H{z) 

for a spatially flat universe. Here, Prob(/Cpeak,oba| Allpeak,h) is given by Eq. (171. In the following, we adopt the model of halo 
mass function in Bhattacharya et al. ]( [MI| ). 


Convergence peak auto spectrum and cross spectrum 


We next consider the auto-correlation function of peaks, and the peak-convergence cross correlation. Marian et al. (20131 
study cosmological information obtained from the statistics using a large set of numerical simulations. They conclude that 
using the auto- and cross-correlation functions can improve the constraints on cosmological parameters when combined with 
the number of peaks. We develop an analytic halo model in order to predict the correlation function of peaks and cross 
correlation between peaks and convergence. 

In the halo model, the number density field of weak leasing selected clusters is given by 
nci(a:) = y^ S^^\x - Xi)S{zi,Mi) 

i 

= J AM S{z, M) Sd{M — Mi) j A^x 5^^\x — Xi)5''^\x — x), (28) 

where S{z, M) represents the selection function (e.g., Takada fc Bridle|2007 1. For clusters identified as leasing peaks, S{z, M) 
is expressed as 


POO 

Si^Z^ M\Uth.re') — I d/Cpeak,obs Prob(A 3 peak,obs | ^peak,/i (-2^5'^^)) • 

" ^thre^noise.O 

Also, underlying matter density field can be approximated as 
pm{x) = y^Phjx- Xi\z,M) 

i 

= J AM M 5d{M — Mi) J A^x'S^^\x'— Xi)um{x — x'\z, M), 


(29) 


(30) 


where ph is the density profile of a massive halo given by Eq. (111, and ph{r\z, M) = Mum{r\z, M). 

In order to derive the auto power spectrum of weak lensing selected clusters for a given threshold Vthie, we first consider 
the auto power spectrum of rid. The two point correlation function of rid is given by 


rickcc(^i^i X 2 ) — 


-b 


(nd(a;i)nd(a;i)) - rid 

S^{Mi)S^^\xi - Xi)S^^\x 2 - Xi)) + { y^ S{Mi)S{Mj)S^^\xi - Xi)S^^\x 2 - Xj)) 
' i,3\i¥=3 

I d^yS\M)5D{M-M,)S^^\x3-y)5^^'>{x2-v)S^^\y-x,)} 

{ y] f AM f A^yS{M)5D{M-Mi)5^^\x3-y)S^o\y-x7) 

i,3\i^3 


X j AM' j A^y'S{M)S d{M'-M j)S‘-^\x2-y')S^^\y' -Xj)) 

= J J - y) 

+ I AM^S{M) I A^y6^^\x3-y) J AM'^S{M') J A^y'S<-^\x 2 - y')^HH{y - y'-, M, M') 

= / I + I dM|^S(M) I AM' ^S{M')^,,Px3-X2-,M,M'),{M) 

where ^hh{y — y'M, M') is the two point correlation function of dark matter haloes with mass of M and M'. In the above 















8 M. Shirasaki et al. 


calculation, we use the following relations as 




{ ^ 5d{M - Mi)5D{M'- Mj)5^^\x^-xi)5f{x2-xi)) 
i,i \ i¥=i 


AM’ 

An An 
dMdM 


■^^hh{xi - X2-, M, M'). 


(32) 

(33) 


We further approximate ^hh as bh{M)bh{M')^ss, where bh is the linear halo bias and is the two point correlation function 
of linear matter density field. Then, we finally obtain the following equation: 


nliicc{xi-X2) = J AM S^{M)S'§\xi - X 2 ) + j AM S{M)bh{M) ^ss{xi-X2). 

The Fourier transform of Eq. (341 is the three-dimensional power spectrum of Ud , which is expressed as 

Pcc{k) = J dVe-‘''-^^,,(r). 

The observed surface number density of weak leasing selected clusters (for a given threshold r'thre) is thus given by 


(34) 


(35) 


(36) 


where A^peak is defined by Eq. (27l. Using the Limber approximation, we can derive the angular power spectrum of p as 

j 2 t. \ 2 






Except for the shot noise term in Eq. (42 I, we obtain 

1 dV\"^ 


[I 


Ht) 

AM —{z, M)S{z, M|i.thre)bh(2, M) 


Pm ( k , Z{x) J , 


(37) 


(38) 


where is the linear matter power spectrum. Throughput this paper, we adopt the functional form of bh proposed in 
[Bhattacharya et al.| ( |2011[ ). 

The similar derivation can be applied to the cross power spectrum of weak leasing selected clusters and leasing convergence. 
Let us first consider the three-dimensional cross correlation function of Ud and Pm'- 


Pm,^cs{xi-X2) = {nd{xi)pm{xi)) - UdPn 

An 


AM — S{M)Mum{xi - X 2 \M) + 


AM^S{M)bh{M) 


Pm^Ssi^i^ 3 ^ 2 )- 


Through the derivation of Eq. (391, we use the following fact that 

pm = j AM ^^hh{M)M J A^x' Um{x - x'\M). 

Then, we can obtain the angular cross power spectrum of p and leasing convergence k by the similar calculation as Eq. 
The cross power spectrum Pp^ is given by 


(39) 


(40) 


S) = / 


W4x) 


1 dV \ 


(x)^ VfVpeakdxdny 




where Pcs{k) represents the three-dimensional cross power spectrum of rid and matter overdensity field 5, i.e.. 


Pcs{k) = J 


= / dVe “^'^Crf(r)- 


Finally, the cross power spectrum between peaks and convergence is given by (also, see Oguri & Takada (20111) 


Pp.w = Pp!:ii) + Pp!:{£), 

J ^ r{xr ViVpeakdxdfly'y 


Pp!:(e) = I 


dx 


r(x )2 V^peakdxdfl 

W4x) /I dV \ 


dM^5(^,M|l/thre) 


M 

pm{z) 


Pm I k — 


r{x) 


z(x),M , 


r(x)^ V-^peakdxdn 
where Um is the Fourier transform of Um{r\z, M). 


) /dM^( 2 ,M)S(«,M|ir,h,e) 6 i.( 2 ,M) 


Pm\k = 


r{x) 


, 2 (x) 


(41) 

(42) 

(43) 

(44) 

(45) 
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2.4 Covariances between statistics 


We summarize covariance matrices between statistics of interest in the following. As a first approximation, we can use the 
Gaussian covariance between three binned spectra Pkk, Ppn and Ppp as, 

1 


Cov[PxY{e),PABie')] = 


{21 + l)At/sky 


[PxmP?B{t) + PxmP^lP)] 5u', 


(46) 


where At is the width of binning in multipole and /sky represents the observed sky fraction. The observed spectra P^yP) 
are then defined by 

_2 . 


pODS _ p I 

KK - K.K “T 


7 pobs _ p j_ 

2ngai ’ iV, 


peak 


pobs _ p 

? -^pK — J pft • 


(47) 


The non-linear gravitational growth causes mode-coupling of the density fluctuations with different wavelengths. The mode¬ 
coupling then induces the correlation of weak leasing statistics between different multipoles, i.e., we can not use the Gaussian 
approximation to covariances (e.g., Sato_et_aJj2009|). Modeling of the non-Gaussian covariances is still being developed (e.g.. 


Gooray & Hu 2001 Takada & Jain 


2004 


Takada & Bridle 20071. We here present a theoretical model of non-Gaussian 


covariance matrices between Ppre and A^peak- 

Let us consider the following set of four-point correlation functions in Fourier space: 

{5m{kl)Sm{k2)5m{k3)Sm{k4)} = (27r)^(5^'(fcl234 )Ti«5« (fel, ^2 , fcs , ^ 4 ) 

(Jc/(fcl)5m(fc2)Jci(fc3)5m(fc4)) = (27r)®(5^’(fci234)T’c«c«(^1, fe2, ^3, fc4) 

{6cl{kl)6m{k2)5m{k3)Sm(k4)) = {2n)^ S^^\kl234)TcSSs{kl, k2, ks, fe4). 


(48) 

(49) 

(50) 


where kij...„ = -f fcj -f ■ • ■ -f fc„, J™ and 5ci represent over density of matter and weak lensing selected clusters, respectively. 
In the flat sky approximation, we can relate three tri-spectra {Tssss,TcScS, Tcsss) with the non-Gaussian part of the covariance 
matrix of weak lensing statistics as follows: 


Gov[p.4^),p.4/)]^^ = ^ I 

Cov[Pp4£),Pp4/)]^^ = ^ I -4,), 

CoY[Pp^ii),PUe')]^o = 

where (f) is the angle between two vectors £ and £'. In practice, the integral over (j) is often simplified as, e.g., 

and are given by 

, IF." 

Jo 


and so on. With Limber approximation, T...., Pp.p., 


c(-£i,-^2,-^3,-^4) = 


dx 


’’(x)' 

IF^ 


;Ts5S5{kl, k 2 , fes, ^ 4 ; z(x))i 


Tp^pPil,£2,i3,£4) = I dx TcSc5{kl,k2,k3,k4-,z{x)), 

Jo r{xr V-^peakdxdH/ 


TpKKK{£l,£2,-t-3,-(-4) — 


r‘ , iF.^ ( 

Jo ^ rixP ( 
Jo ^ rixP ( 


dV 

r(x)® V-^peakdxdfl 


Psss{ki,k 2 , k3, k4-z{x)), 


(51) 

(52) 

(53) 

(54) 

(55) 

(56) 

(57) 


where ki = L/x> Xs is the comoving distance to sources, the window function IF. is given by Eq. (261 and Apeak is defined 
by Eq. 

Previous works show that the dominant contribution of the non-Gaussian covariance is the so-called one-halo term of the 
relevant tri-spectrum at £ 100 (e.g., Sato et al.|2009 1. One halo term arises from the four point correlation between different 

modes ki {i = 1, 2,3,4) in a single dark matter halo. Thus, one halo term of the underlying tri-spectra between Sm and Sc can 
be calculated as 

4 

u„(fci|2;, M)Umik2\z, M)Um{k3\z, M)Um{k 4 \z, M), 


^5i5i5i5(^ 1 5 ^25 ^37 ^4 7 ■S) 

/ ,, . dn 
j AM 

^c5ci5(^1 7 ^2 7 ^37 ^47 ■^) 

j AM 

^ci5i5i5(^ 1 7 ^27 ^37 ^47 ■^) 

J AM 


{z,M) 

{z,M) 

{z,M) 


M 


pm {z^ 

M 

pm {z^ 

M 

pm {z^ 


S{z, M)Um{k2\z, M)Um{k4\z, M), 


S{z, M)Um{k2\z, M)Um{k3\z, M)um{k4\z, M). 


(58) 

(59) 

(60) 


Another important contribution of covariance at degree scales or less is so-called halo sampling variance (HSV) (Sato et al. 
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2009 Kayo, Takada fc Jain|[20l3l. This term describes the mode-coupling between the measured Fourier modes and larger 


modes of length-scales comparable to survey volume. It is expected to be important when the number of massive haloes found 
in a hnite region is correlated with the overall mass density fluctuation in the region (Hu fc Kravtsov|2003 1. Following Kayo, 


Takada & Jain (20131, we model the HSV of weak lensing statistics as 

Cov[P.«W,;’.4^')]hsv = [|dM^6,(M)|S,MM)r 




f 

JO 


hc\k 

^P,^(fc)|tK(fcXesu.vey)|" 


C0v[Pp4f),Pp4/)] 


Cov[Pp„(^),p«4/)] 


HSV 


HSV 




^P,^(fe)|lK(fcX©survey)|" 


AM —hH{M)S{z, M)k4£|M) 


AM' ^hH{M')\kH{t'\M't 


f 

JO 


^P,^(fc)|tK(fcX©survey)|' 


(61) 


(62) 


(63) 


where kh is the Fourier transforming of Eq. (131. Here, W (fcxQaurvey) represents the window function of the survey region in 
Fourier space with ©survey denoting the squared root of the survey area. We use the circular function of W{x) = J\{x)/x. 

Cross covariance between the number count of weak lensing selected clusters and the lensing spectra can be naturally 
incorporated in the halo model ( Takada k, Bridle|2007 Takada fc Spergel|2014 1 . 

The covariance of A^peak with two different thresholds is given by 

Cov[Arpeak(r'thre,l),fVpeak(!^thre,2)] = f ^ ^ f ) f f 

X [riclixO, z(x)WthYe,l)ncl{x , z{x')Wthre,2) - hcl(2:(x)kthre.l)hcl(2(x0kthre,2)] ,(64) 


where W{6) represents the window function in real space. Performing the similar calculation as in Eq. (311, one can find that 
(see also the appendix in Takada & Bridle ( 2007| >) 


Cov[A^peak(^thre,l)7 -^peak (i^thre,2)] — *^12 


-^p eak ( r't hre, 1) 


47r/, 


sky 


/ 


+ / dx 


VdxdHy' 


r{x) ^ 

f 


27r 


/ dA^^6h(A^)5'(2,Mkthre,l) 

Ph^mw{kxQ.^Y..y)\ 


(65) 


where /sky is the sky fraction for survey of interest. In order to derive the cross covariance between A^peak and P„„ or Pp^, we 
hrst define the estimator of power spectrum as 




( 66 ) 


where the summation is taken over all the Fourier modes in the range of [f —At'/2, £-|-Af/2] and A^ is the width of multipoles. 
Also, Np{£) represents the number of modes to estimate of power spectrum with the multipole of £. With Eq. ( | 66 [ ), we can 
express the cross covariance of Apeak and P^k as 


Cov[Apeak(r'thre), P,^,^(^)] = 


4nfskyNp{£) 


^ / A^eW{e){R{£)R{-£) [nci(x0,z(x)kthre) - nci(2(x)kthre)]), 


47r/skyAp(£) 


A^e tK(6») Apeak(l'thre) ^ ^ (k(^)k(-^) p(£>thre))6“' 
£ £' 


£'■0 


(67) 


where we use Eq. (36 1 through the derivation. Similarly, the cross covariances of Apeak and Pp^ is given by 

Cov[Apeak(r'thre,l), PpK(^kthre,2)] = “j- j / d^S W(0)Apeak(l^thre,2) 

47r/skyAp(tj J 


X T T(P(^l^thre,l)K(--g)p(^Vthre,2))e 
£ £' 


-i£'-6 


( 68 ) 


Therefore, the cross covariances of Apeak and lensing spectra include the three point correlation of the relevant field p or k. 
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Lbox [A-iMpc] 

•^init 

No. of sim. 

output redshift 

450 

72 

10 

0.025, 0.076, 0.129 

900 

36 

10 

0.182, 0.237, 0.294 

1350 

24 

10 

0.352, 0.412, 0.475 

1800 

18 

10 

0.540, 0.607, 0.677 

2250 

15 

10 

0.751, 0.827, 0.901 

2700 

12 

10 

0.990, 1.077, 1.169 


Table 1. Parameters used for A^-body simulations. Each simulation was run with 1024^ dark matter particles. The output redshift 
of each simulation corresponds to the comoving distance to the center of lens-shells. We adopted the standard ACDM model, which is 
consistent with WMAP nine-year results ( [Hinshaw et al.|2013| l. 


We can thus summarize the relevant covariance are 

C0v[iVpeak(;^thre),P«4^)l = ^ ^ (0, ^/X, ^/X; ^(x)), (69) 

Cov[Arpeak(^'thre,l),r’p«(^kthre,2)l = j ■B^{l)c(2)« (0, ^/X, ^/X; ^(x)), (70) 

where Bds and -Bc(i)c( 2)5 are defined by 

(5cl{kl)S^{k2)5m{k3)) = (27r)®(5®(fci23)-BcM(fcl,fc2,fe3), (71) 

(Eel (^11 i^thre,! (^2 | ^^thre,2)^ m (fc3)) = (27r)®5®(fel23)B c(l)c(2)5(fcl, fc2, fcs)- (72) 


At I ^ 500, the main contributor to the covariances is the one-halo term of the above bi-spectra. Similarly to the case of 
tri-spectra, the corresponding terms are expressed as 


Blss{ki,k2,k3-,z) 


Bl(i)c( 2 )s{ki,k 2 , fc3; z) 




M 

pm (^) 

M 

Pm{z) 


S{Z, M|r'thre)Mm(fe2|2, M)Um(k3\z, M), 
S{Z, M\UthTe,l)S{z, M\Vthie,2)Um{k3\z, M). 


(73) 

(74) 


3 NUMERICAL SIMULATION 

In order to study in detail the weak lensing statistics considered in this paper, we use mock weak lensing catalogs generated 
from a set of full-sky weak gravitational simulations. The A-body simulations reproduce the gravity-driven non-Gaussianities 
in the underlying matter density field, and realistic mask regions are pasted on our mock catalogues. Thus we can compare 
our model directly with the simulated catalogues. 


3.1 A-body simulation 


We first run a number of cosmological A-body simulations to generate three-dimensional matter density fields. We use the 
parallel Tree-Particle Mesh code Gadget2 | Springel||2005 1. We arrange the simulation boxes to cover the past light-cone of a 
hypothetical observer with an angular extent of 47r/8 steradian. A set of simulations consist of six different boxes and cover 
one-eighth of the sky. In order to cover the full-sky, we use the same boxes by adopting periodic boundary condition (see 
Figure 2). 

The simulations are run with 1024^ dark matter particles in six different volumes: the box side length ranges from 
450/i“^Mpc to 2700/i“^Mpc with increments of 450/i“^Mpc. The largest volume simulations with 2700/i“^Mpc on a side 
enable us to simulate the gravitational lensing effect with source redshift of ~ 1. We generate the initial conditions using 
a parallel code developed by Nishimichi et al. (20091 and Valageas fc Nishimichi] (20111, which employs the second-order 
Lagrangian perturbation theory (e.g. Crocce, Pueblas fc Scoccimarro||2006 1 . We set slightly different initial redshift 2init as 
the box size increases. In order to generate the initial conditions, we calculate the linear matter transfer function using GAME 


(Lewis, Challinor & Lasenby 20001. Our fiducial cosmological model is characterized by the following parameters: matter 


density flmo = 0.279, dark energy density Hao = 0.721, the density fluctuation amplitude ag = 0.823, the parameter of the 
equation of state of dark energy wo = —1, Hubble parameter h = 0.700 and the scalar spectral index ris = 0.972. These 
parameters are consistent with the WMAP nine-year results (Hinshaw et al. The parameter of our A-body simulations 

are summarized in Table [T] 
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Comoving radial distance [Mpc/h] 0.129 0.294 0.475 

Redshift 

Figure 2. The configuration of spherical shells of projected matter density. The left panel illustrates how we place the outputs of the 
Ai-body simulations. There, the dashed line represents the original box size and each of the three shells is taken from the snapshot at 
different redshift. The right panel shows the nested structure of simulation boxes and the configuration of multiple lens shells. In the 
right panel, different colors are used to indicate different sets of simulations. The corresponding redshift of each shell is also shown at 
the bottom of the right panel. 


3.2 Ray-tracing simulation 


We briefly summarize our ray-tracing simulations with full-sky coverage. The detailed description of our ray-tracing simula¬ 
tions is found in Appendix]^ In our ray-tracing simulation, the light ray path and magnification matrix are calculated by 
the standard multiple lens-plane algorithm. The multiple lens-plane algorithm on a spherical geometry requires contiguous 
spherical shells of projected matter density. We thus utilize A^-body simulations in Section [3.1| to generate a number of thin 
shells with width of 150/i“^Mpc and produce three shells from a single simulation box as shown in the left panel of Figure 
In order to extract a target shell region on a lightcone, we choose an appropriate snapshot at the redshift that correspond 
to the comoving distance to the shell from a hypothetical observer point (origin). A set of projected density shells are then 
configured using the nested structure of the simulation boxes. The right panel of Figure shows the configuration of shells 
used in our ray-tracing simulations (only for a 0.5). 

In this paper, we use HEALPix libraries for pixelization of a sphere | |G6rski et al.|200^ . The angular resolution parameter 
nside is set to be 4096. The corresponding angular resolution is ~ 0.86 arcmin. For each shell, we calculate the projected mass 
density map from A-body particles using Nearest Grid Point method. Then, we perform the spherical harmonic transformation 


of the density shells to calculate the gravitational lensing potential on each shell via the Poisson equation (see, Eq. (C4|). The 
obtained gravitational potential and its derivatives are used in the standard multiple lens-plane algorithm. In our simulation, 
we follow light ray trajectories from z = 0 to z = 1. The position of a ray and the magnification matrix Aij at each shell are 


updated by the recurrence relation as in e.g., Hilbert et al. (20091. The initial position of each ray is set to be the center of a 


pixel on HEALPix map. As a ray propagate with deflection, the position at each shell deviates from the position of the HEALPix 
map in general. We evaluate the potential at each shell by using an inverse-distant weighted interpolation of the potential 


field. The formalism of the multiple lens-plane method on a sphere is found in e.g., Teyssier et al. (20091; Becker (20131. 


We obtain a total of ten all-sky convergence maps from ten sets of realization of A-body simulations. Note that we use 
different initial random seeds to generate A-body simulation data in order to avoid finding similar structures along a line of 
sight through ray-tracing. Figure]^ shows an example of our simulated convergence map on a full-sky. 


3.3 Shape noise and sky mask 


We generate realistic mock catalogues by adding the effect of shape noise and sky mask to the obtained convergence maps. 
The intrinsic ellipticity of source galaxies is the main contaminant of cosmic shear measurement. We model the shape noise 
in a convergence map by assuming a two-dimensional Gaussian field as follows: 


{K!Ss{di)liN{dj)) = 


2rigalflpix 


( 75 ) 


where a-y = 0.4, rigai = 10 arcmin and flpix = 47r/12/4096^ = 6.24 x 10 ® str. 
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Figure 3. One of our simulated full-sky convergence maps. We perform the Gaussian smoothing with smoothing length of 2.12 arcmin. 
The shape noise is not included in this figure. The red regions represent high convergence, while the blue regions correspond to negative 
convergence. 



RA [deg] 



Figure 4. The configuration of the masked regions. In the left panel, black regions show the masks on the bright stars with the R-band 
magnitude of Mji < 10. The right panel represents the typical configuration of the masks in a circular region with radius of 60 arcmin. 
Black isolated regions show masked regions due to the bright stars and white regions correspond to additional masked region. These 
white regions remove some pixels affected by the convolution of survey masks (black regions) with the smoothing filter. 
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It is important to use a realistic sky mask to study statistics for upcoming lensing surveys. We here consider masks 
owing to bright stars. Among the planned HSC survey regions, we select two continuous regions with sky coverage of ~ 565 
and ~ 680 squared degrees. J2000 coordinate of these regions are given by 22'’00™ < RA < 2''40™, —1° < DEC < +7° and 
8’'30™ < RA < 15’'00™,—2° < DEC < +5°. In the following, we consider the bright stars with the R-band magnitude of 
Mij < 17. We then select 1,149,871 stars located in these two regions from USNO-A2.0 catalog In this paper, we assume 
the following relation between the R-band magnitude Mr and the effective radius of a halo of bright star Tatar: 


Tatar [arcsec] = 


180 
0.2 X 


{Mr < 9), 
{9<Mr< 17). 


(76) 


Then, we paste a circular mask with Tatar around the pixel located at each star using Eq. (76l. If Tatar is less than the angular 
size of our all-sky map (~ 1 arcmin), we simply mask the pixel located in the star. After the above procedure, we further 
remove the “isolated” pixels whose surrounding pixels are all labeled as masked pixels. The final mask configuration generated 
in this way is shown in Figure The left panel displays the masked region over the 565 -I- 680 square degrees. The black 
dots shows the masked pixels on very bright stars {Mr < 10) and smaller masked regions are distributed in the white region 
homogeneously (but not shown in the figure). The mask covers over ~ 280 squared degrees in total. The right panel shows 
an example of our masked sky simulation in a circular region with a radius of 1°. 

From a single full-sky simulation, we make 20 mock HSC convergence maps with the masks by choosing the desired sky 
coverage (565 -I- 680 = 1245 squared degrees). We allow to have small overlap regions between the 20 masked maps. The 
overlap regions are located near the edge of the HSC sky coverage, and thus we expect this minor compromise does not affect 
the final results significantly. Finally, using 10 independent full-sky maps, we obtain a total of 20 x 10 = 200 realizations of 
mock HSC lensing catalogs. 

For a smoothed convergence map, we apply the different masking. As shown in Eq. ( |76[ ), stars with Mr > 12 have smaller 
mask radii than the pixel size in our simulation. Therefore, in principle, we can extract some information from the pixels 
where stars with Mr > 12 are located. Because such faint stars would not affect the smoothed convergence map, we mask 
only the bright stars with the R-band magnitude of Mr < 12. There are also ill-defined pixels that are compromised by 
the convolution with the smoothing filter. We remove such ill-defined pixels within 5 arcmin from the boundary of the mask 
regions. Consequently, masked regions on the smoothed convergence tC maps differ from the original masked regions. Finally, 
we have a total of ~ 412 squared degrees as unmasked regions. The corresponding sky fraction is ~0.01. 


4 RESULT 

We present the weak lensing statistics measured from the full-sky simulations and also those from our 200 mock HSC 
catalogues. In the following, we define the threshold of lensing peak as K./a^a\se,o, where (Tnoise.o is given by Eq. ( |10[ |. Note 
that we use this definition also in the case without noise. 


4.1 Ensemble average of statistics 

All-sky 

We first show the ensemble average statistics over the ten full-sky simulations. They can be regarded as the expected values 
from an idealized full-sky observation. Figure summarizes the measurements of Pkk, Ppp, PpK and A^peak. We selected the 
lensing peaks with threshold of vthre = 3 (or tC ~ 0.05) in both the maps with and without shape noise. To calculate the 
correlation in harmonic space, we correct the pixelisation effect with the pixel window function of HEALPix. Overall, our model 
is in good agreement with the measurement from the full-sky simulations. In the case without shape noise, we calculate Ppp, 
Ppn and Apeak assuming the one-to-one correspondence between lensing peaks and dark matter halos, i.e., 

Prob(/Cpeak,obs| Apeak,/i(^i Af)) ~ (Apeak,obs Apeak,(^; Af)) . (^'^) 

The results of Ppp and Pp^ from the maps with and without noise show appreciable differences at large angular scales (£ 100). 

This can be explained by the modulation of peak height due to the shape noise. When we select the lensing peaks with a 
given Vthie, we effectively include less massive dark matter haloes as well in the noisy A maps. Then both Ppp and Pp^ are 
reduced at large angular scales in the noisy maps because less massive halos have weaker clustering. Also, the number count 
of peaks in the noisy convergence maps is described well by our model as shown in Section [2.3| In the maps with shape noise, 
we detect more lensing peaks for a given threshold, for example, by a factor of ~ 3 for fthie = 3. 


4 


http: //tdc- WWW. harvard.edu/catalogs/ua2.html 
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Figure 5. We compare the statistics measured from the ten all-sky simulations and our model prediction. In each panel, the black line 
is our halo model prediction. The red points with error bar represent the measured signal from the ’clean’ convergence maps without 
noise. The blue points show the result with noise. In the bottom portion in top right and the bottom left panel, we show the ratio of the 
two. The error bars indicate the standard deviation over ten realization. Note that the threshold of lensing peaks is defined by /C/crnoise,0 
regardless of the presence or absence of noise. 


Masked sky 


In the presence of masks, we need the correction for the mode-coupling effects in measurements of power spectrum. We adopt 


the pseudo-spectrum estimator for this purpose (Hansen & Gorski 2003 Efstathiou 2004 Brown, Castro &: Taylor 2005 


[Hikage et al.|20rT| ). The basics of the method is summarized in Appendix D In measurements of the binned two spectra 


we follow the similar method shown in Brown, Castro Sz Taylor (20051; Hikage et al. (20111. Let us consider the 5-th 

min < ^ < ^min- casc, we can calculate the band-powers (Pb) as 


and Pj 

binned power spectrum in the multipole range of 

b' i 

where P(^) is the measured power spectrum on a masked sky and tSbe is the binning operator, which is defined by 

£{£ + l)/(2^)/(4+l - 4i„) (4in < t < 4tn), 


Bit = 


0 (otherwise). 

Here, we define the binned coupling matrix Mai as 
Mbb' = BbiMui Qijiy , 


(78) 


(79) 


( 80 ) 
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Figure 6. The auto- and cross-power spectra obtained from 200 mock HSC catalogues. The black line is the result of our halo model. 
In the left panel, the dashed line shows the shape noise contribution to convergence power spectrum. 


where the definition of M^/ is found in Appendix [D] and is given by 

/ 2VW^+l)l (Cn<^<4tn), 

I 0 (otherwise). 


(81) 


The number of bins is set to be 30. We perform the binning in linear spacing for the first 10 bins as = 2 + 106 (1 < 6 < 10), 
while the remaining 20 bins have logarithmically eqnal spacing up to £niax = 2000. The measured and corrected binned power 
spectra are plotted in Figure We use 200 masked sky simulations to calculate the average of the binned and Ppn- 

For cross-correlation, we find ~ 100 — 200 leasing peaks with t'thre = 3 on each simulated HSC map. In Figure]^ the blue 
points with error bars shows the measured spectrum, and the red points are the corrected spectrum with the pseudo-spectrum 
method. Clearly, we can recover the underlying power spectrum with the correction of the mode-coupling due to masks. The 
difference in amplitude can be explained approximately by the effective fraction of sky coverage (~ 0.01 — 0.02) 


4.2 Covariance 


We use 200 masked sky simulations to calculate covariances between the weak leasing statistics of interest. First, we calculate 
the convergence power spectrum covariance. In the flat sky approximation and without masks, the covariance of the binned 
power spectrum is expressed as (e.g., Cooray fc Hu|20'oT I 


Cov[p^4ei),p^4£,)] = 


47r/sky Ji. As,i Jn As,j 


(82) 


where Pnniit) represents the i-th binned power spectrum, /sky is the fraction of sky covered by the observation, A£ is the bin 
width in £ space, and As,i is the area of the two-dimensional shell around the i-th bin £i in Fourier space. Here, T^kkk is the 
tri-spectrum of convergence, of which definition and modeling are found in Section |2.4| In practice, we simplify the second 
term in Eq. (|82|) as 

(83) 


/ -I,£',-£') 

J Jl- AsJ 


In the presence of masked regions, the estimator of the binned power spectrum is given by more complex expressions given in 
Eq. ( |78[ ). The mode-coupling due to masked regions induces the intricate correlation between different Fourier modes in the 
covariance of Pkk (the exact expression is found in Appendix 0- Nevertheless, in the case where the value of masked pixels 
is set to be either 0 or 1, we can use the simplified formula of covariances as shown in Efstathiou (20041. When we work with 
harmonic space, the approximated formula can be expressed as 


Cov[P.44),P«4^2)]„.ask E 


f KK 
p! p! 


2P^4£'4P^4£'4 + — V(2Pi + 1)(2£[ + 2)X««4/i,/i,4,£^) 

47r 


where represents the mode-coupling matrix for convergence. In Eq. (841, we evaluate the term Tkkkk by the tri-spectrum 
in harmonic space following Hu (20011 because we define Tkkkk in the fiat sky approximation. 


, (84) 
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Figure 7. The variance of the psudo-spectrum estimator of Pk.k measured from 200 masked sky simulations. The black points show 
our measurement. The three lines represent the different contribution to variance. The black line is the simplest halo model with the sky 
fraction of /gjjy = 0.023, while the red lines show the halo model covariance with the correction of the effect of masked regions. The red 
dashed line corresponds to the Gaussian model including the effect of masked region. For the red solid line, we take into account both 
non-Gaussianities caused by gravity and masked regions. 


We can then examine the validity of Eq. (821 by using the measured covariance of Pk,k over 200 masked sky simulations. 


We use the same binning as in Section [4.1| but reduce the number of bins to 10 by taking average of the binned powers over 
nearest I bins. Figurej^shows the measured variance of the pseudo-spectrum estimator of Pkk - The black point shows the result 
obtained from the 200 simulations and the solid line represents our model of covariance in Eq. (82 I with the appropriate value 
of /sky for our simulations. Interestingly, although the simple model of Eq. (821 is expected to account for non-Gaussianities 
caused by gravity, it underestimates the actual covariance by a factor of ~ 10. The corrected covariance components are shown 
in Figure]^ The first term and second term in Eq. (841 are plotted as red dashed and red solid line, respectively. The overall 
amplitude of the variance can be explained by the mask correction, while the contribution from tri-spectrum dominates at 
£ 200. Clearly, it is problematic to adopt the commonly used estimate of covariance given by in Eq. (82l. 

To perform combined analysis of Pk,k.,Ppk and A^peak, we need cross covariance between the statistics. Unfortunately, it 
is difficult to model the impact of masked region on the cross covariances in a similar manner to Eq. (84|. Here, we simply 
compare the measured cross covariances with the prediction of our halo model. In the model, the covariance can be derived 
similarly to Eq. (82l. For example. 

Si, 
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where we have defined Tp^pK, Tp^KK, and the other covariances (i.e. Cov[Apeak, A^peak], Cov[Apeak, Pkk], and Cov[Apeak, Pp«;]) 
in Section f2.41 

Figure]^ shows the measured peak-power cross-covariance. In constrast to the case of Pk,k, the covariance including Pp^ 
and Apeak is less affected by masks: the simple model that accounts for /aky yields a reasonable result with respect to that of 
the simulations. We find that the difference is by a factor of ~ 3 at most. 


5 CONCLUSION AND DISCUSSION 

We have performed all-sky lensing simulations to generate a large set of realistic weak lensing mass maps with complex masked 
regions by incorporating the actual position of bright stars. We have used the set of the mock samples to study the statistical 
properties of weak lensing convergence and convergence peaks in detail. Full nonlinear covariances between the statistics have 
been also calculated from 200 realization of masked maps. We have also developed an analytic halo model that provides 
reasonably accurate prediction for the statistics. 

When adopting a Gaussian smoothing with the full width at half maximum of 5 arcmin, we can associate weak lensing 
convergence peaks with dark matter halos with mass of ~ 10^"^ h~^M q at z ^ 0.1 — 0.2. We can also estimate the modulation 
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Figure 8. The cross covariance of weak lensing statistics calculated from the 200 simulations. In each panel, the black or colored points 
show the measured covariance, and the solid line shows the corresponding halo model prediction. We define the peaks with r'thre when 
measuring PpK- In the bottom two panels, the color indicate the results for different threshold r'thre in ^peak- The threshold value is 
varied in the range t'thre = 3.5 — 7.7. 


of peak height due to shape noise by using a model based on Gaussian peak statistics. Thus, the abundance of peaks, the 
angular correlation function, and the cross-correlation of peaks and cosmic shear are all obtained accurately by our halo model 
approach with the shape noise correction. Furthermore, the halo model can also take the mask effect into account and indeed 
produces accurate ensemble average of statistics and their covariances. The impact of masked regions on the covariance can 
be described by the following two effects: (i) reduction of sky coverage and (ii) the mode-coupling effect between different 
Fourier modes. We find that the former affects the overall amplitude of cross-covariance between statistics, while the latter 
is important for the covariance of cosmic shear power spectrum Pkk- For the masked regions adopted here, ignoring the 
mode-coupling effect would induce underestimation of the covariance of P^k, by a factor of ~ 10(!). 

The number density of source galaxies is an important factor in the statistical analysis of weak gravitational lensing. As 
one may expect, a large number density of sources is desired to find clusters with high accuracy and perform cosmological 
analysis with selected clusters. In comparison with the case of rigai = 10arcmin“^, we have confirmed that the signal-to-noise 
ratio increase by a factor of ~ 1.5 in combined analysis with Pk,k., Ppn and A^peak even if we ignore the shape noise contaminant 
(i.e., rigai —>■ oo). This suggests that imaging over a wide area is suitable for cosmological analysis with the lensing statistics 
even if the number density of sources is not significantly increased in such surveys. 

We have also examined the validity of our model for two additional cases: Ugai = 5 and 30 arcmin~^. When the smoothing 
scale, the rms of intrinsic ellipticity, and the source redshift are all hxed, the halo model prediction is in good agreement with 
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the result of our full-sky simulations in the case of rigai = 30arcmin“^, but the agreement is worse with ngai = 5arcmin“^|^ 
Therefore, when we consider the typical value of the rms of intrinsic ellipticity and the source redshift, our model is expected 
to be accurate when rigai ~ 10arcmin“^ with ~ 2 arcmin Gaussian smoothing. 

Our model has been successfully applied the statistics of a smoothed convergence map with shape noise. In principle, one 
can hnd an optimal filter function so that the number of detected clusters is increased (e.g. Hennawi &: Spergel|2005 Maturi 


|et al.|[^005[ ). Our model is based on the assumption that shape noise in a smoothed lensing map has Gaussian properties. 
The assumption is valid when the ellipticities of the source galaxies are uncorrelated and when there are a sufficient number 
of source galaxies per pixel (i.e. the central limit theorem). We expect that, while our model can be applied to the general 
form of Hlter function, it may not be appropriate when there are only few source galaxies and/or when there is no one-to-one 
correspondence between peaks and halos. 

Clusters of galaxies are important targets in future cosmological surveys. Selecting galaxy clusters based on cosmic 
shear measurement does not rely on mass estimate nor on calibration with additional information such as X-ray brightness. 
Nevertheless, there remain some systematic effects worth mentioning here. 

Gravitational lensing causes not only distortion of source images but also magnification. Using magnified galaxies in a 
flnx- and size-limited survey would potentially cause systematic effect (s) in the weak lensing statistics. The magnification 


effects on lensing peak statistics have been already studied in e.g., Schmidt & Rozo (20111. Recent numerical study by Liu 


|et al.| ( [2014a[ ) suggests that the magnification effect causes non-negligible bias in parameter estimation in the case of LSST. In 
order to examine the magnihcation effect further, it is essential to run high angular resolution simulations. This is because the 
mean number of source galaxies on each pixel should be less than unity to make the one-to-one correspondence between a pixel 
and a (magnified) source galaxy. In the case of rigai = 10arcmin“^, we should set the pixel size to be 1 /.^rigai ~ 0.3 arcmin. 
We will perform such simulations to study the magnification effect in wide-held surveys in detail. 

Another important issues are uncertainties and systematic bias associated with baryonic effects. Previous studies (e.g.. 


Semboloni et al. 2011 Semboloni, Hoekstra & Schaye 2013 Zentner et al. 20131 explored the impact of the baryonic component 


to two-point statistics of cosmic shear and consequently to cosmological parameter estimation. The baryonic effect is likely 
important in weak lensing peak statistics. Indeed, Yang et al. (2013) show appreciable baryonic effects on peak statistics using 
a simple model applied to dark-matter-only simulations, whereas the baryonic effect on higher order convergence statistics 
have been studied with numerical simulations ( |Osato, Shirasaki fc Yoshida|2015| | . Recently, [Mohammed et aT)] ( |2014| ) explored 
halo model approach to include the baryonic effect on cosmic shear statistics. 

The statistical properties and the intrinsic correlation of source galaxies and lensing structures are still uncertain but 


could be critical when making a large lensing mass maps. Among such correlations, source-lens clustering (e.g., Hamana et al. 


2002) and the intrinsic alignment (e.g., Hirata fc Seljak|2004) are likely to compromise cosmological parameter estimation. A 


promising approach in theoretical stndies would be associating the source positions with their host dark matter halos on the 
light cone. This is along the line of our ongoing study using a large set of cosmological simulations. 

Weak gravitational lensing is a promising tool to probe the dark matter distribution in the universe. Statistical analysis of 
a reconstrncted mass map can be performed to extract precise cosmological information. The peak statistics considered in the 
present paper contain the information related to massive objects such as clusters of galaxies and thus have a great potential 
to probe cosmology and constrain the model of structure formation simultaneously. Ongoing/upcoming imaging surveys such 
as HSC, DES, and LSST in the near future, will provide the largest dark matter map we have never seen before. We expect 
our study presented here provides a useful guide to interprete properly the reconstructed mass map and to reveal the nature 
of the dark components in the universe. 
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APPENDIX A: HALO-PEAK MATCHING 

In this appendix, we examine the correspondence between dark matter halos and the local maximum in lensing convergence 
map. 


® We expect that the disagreement for small rigai would be caused by the offset between the position of a peak and the center of the 
corresponding halo. The offset effect would be more important when shape noise increases as shown in Fan, Shan & Liu 12010 I. 
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Figure Al. The correspondence between dark matter halos and lensing peaks. The top panel shows the scatter plot of peak height and 
the expected convergence by the matched halos in absence of noise. The lower panel corresponds to the case with shape noise. In each 
panel, the horizontal axis represents the peak height and the vertical axis shows the expected convergence of NFW halos. 


In order to generate mock halo catalogs, we identify dark matter halos in outputs of our A-body simulation (see, Sec¬ 
tion 3.11 using the standard friends-of-friends algorithm with a linking parameter of 6 = 0.2 in units of the mean particle 
separation. We use dark matter halos with mass greater than 10^^ h~^MQ in the following analysis. This lowest mass corre¬ 
sponds to the mass of 20 particles in the largest simulation. Then, the position of dark matter halo in A^-body simulations 
are arranged in the same way as in the ray-tracing experiments in Section [3.2[ 

With our ray-tracing simulations and mock halo catalogs, we study the correspondence between halos and the peaks in 
weak lensing convergence maps. We first identify the local maxima in the smoothed lensing convergence field with source 
redshift of ^source = 1. In this appendix, we again adopt the Gaussian smoothing with the full width at half maximum of 5 
arcmin. When including the shape noise in convergence maps, we set — 0.4 and rigai = 10arcmin”^. For selection of peaks, 
the threshold of peak height is set to be /C = 0.03. This value corresponds to ~ 3cr in smoothed convergence maps without 
noise. For a given position of lensing peak, we search for the matched dark matter halos within a radius of 5 arcmin from 
the peak position. This search radius is set to be larger than the smoothing scale but still smaller than the angular size of 
massive halos at « ~ 0.1 — 0.2 (also see, Hamana, Takada fc Yoshida||2004 1. When we find several halos in search radius, we 
regard the matched halo as the closest halo from the position of peak. For each matched peak, we estimate the corresponding 
convergence by using the universal NFW density profile (see Section [2.2| in detail). In the calculation of expected convergence 
from FoF halos, we simply assume that the FoF mass is equal to the virial mass. In total, we find 632,238 and 1,404,538 pairs 
of peaks and halos over 10 noise-less maps and noisy maps, respectively. 

Figure [AT] shows the scatter plot of peak height in K. map and the expected convergence by NFW halos. The horizontal 
axis corresponds to peak height, while the vertical axis shows the corresponding convergence expected by NFW halos. Thus, 
the color map in each panel shows the probability of Eq. 0- We present the line oi y = x as the dashed line in each panel. 
In lower panel of this figure, we show the effect of the modulation of peak height as the solid line with error bars. The solid 
line is derived by ^peak,obs( 2 , A4) in Eq. (22 I and the error bars reflect the scatter of ^peak,obs(-2, ihf). As shown in previous 
works, we confirm the good correspondence between the matched dark matter halos and lensing peaks in the noise-less maps. 
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Figure Bl. The property of lensing peaks on smoothed maps by the compensated Gaussian filter. The left panel shows the scatter 
plot of peak height and the expected convergence by the matched halos in absence of noise. In the left panel, the dashed line shows 
one-to-one correspondence. The purple points represent the mean relation of the expected convergence and measured peak height and 
the best-fit linear relation is expressed as the purple line. The right panel shows the comparison the peak abundance measured from 
ten full-sky simulations and our model prediction. In the right panel, the black line is our halo model prediction. The red points with 
error bar represent the measured signal from the ’clean’ convergence maps without noise. The blue points show the result with noise. 
We normalize the measured peak height by crnoise,0 both cases with and without noise. The error bars indicate the standard deviation 
over ten realization. In the calculation of peak count, we correct the biased relation of /^nfw ^peak shown in the purple line in 
the left panel. The details of correction are found in the text. 


Also, our model as shown in Eq. (221 can explain the average relation between peaks and dark matter halos even in the case 
with noise. 


APPENDIX B: THE CASE OF COMPENSATED GAUSSIAN FILTER 


We examine another practical filter function to construct smoothed lensing maps. In this appendix, we consider Zs, 
and we set cr^ = 0.4 and rigai = 10arcmin“^ when including the shape noise in convergence maps. 

We first consider the compensated filter Uc based on the Gaussian form as 




1 I 

l-exp( 


= 1 . 


(Bl) 


where 9o represents the boundary of the filter and we set Uc to be zero for 9 > 6o- We adopt the smoothing scale of 
9g = 5/y^ 2 arcmin and 9o = 30 arcmin. For the compensated filter function of Uc, the noise power spectrum on a 
smoothed lensing map is expressed as (van Waerbeke|2000| 


pj^p) = 

ingal 


(B2) 


where is the rms of the intrinsic ellipticity of sources, rigai represents the number density of source galaxies, and Uc is the 
Fourier transform of Uc- As shown in Eq. (101, the noise variance Unoise.o is evaluated by the integral of Pm(1) in Fourier 
space. In the case of 9o = 30 arcmin, the boundary of the filter changes Unoise.o by about 1% compared to the case of the 
usual Gaussian filter. Thus, we can safely ignore the difference of Unoise.o between the compensated and the usual Gaussian 
filter for our parameter choice. 

In order to investigate the effect of the modification of filter function on lensing peak statistics, we first study the 
correspondence between peaks and halos by the halo-peak matching analysis as shown in Appendix]^ When we limit lensing 
peaks with the height larger than 0.02, we find 1,045,291 matched pairs over ten full-sky maps without shape noise. The left 
panel in figure [Bl] shows the scatter plot of the measured peak height and the expected convergence signal for the spherical 
NFW halo. In the calculation of expected signal, we simply assume that FoF mass of halos is equal to the virial mass and use 
the model of concentration parameter Cvir in |DufIy et al.| ( [2010[ ). We expect that the one-to-one correspondence between peaks 
and halos would still hold. However, we find biased relation between the mean measured height and expected one even in the 
absenee of noise. The biased mean relation is shown by the purple points in the left panel in figure |B1| It can also be fitted 
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well by linear relation of /Cpeak = o/Cnfw + /3- We find the best-fit value of a is 0.9 while the offset /3 can be approximated to 
be zero. A similar relation is also found by Hamana et al. (2012 I. It might be caused as a consequence of various effects such 
as the mismatch of FoF mass and virial mass dehned by spherical over-density. 

Even without knowing the origin of the biased relation, we can still predict the lensing peak statistics with the compensated 
hlter by adopting the biased relation in our halo model. Our approach is simply to replace /Cpeak,h(2, M) with a/Cpeak,jt(2, M)-\- 
P for the calculation of Eq. (17l. Here, we also assume that ICpea.k,h can be evaluated by 

iCp,^kAz,M)= f Aeu,ie-,OG,eo)Kh{e\z,M), (bs) 


where Kh{d\z, M) represents the convergence profile of spherical NFW halos with mass of M and the redshift of 2 . With the 
above correction, our model can provide a reasonable fit to the measured peak statistics from ten full-sky maps as shown in 
the right panel in figure |B1| The right panel shows the measured peak count for the compensated filter with and without 
shape noise. In the right panel, the black solid (dashed) line represent our halo model in absence (presence) of noise with the 
correction for the biased relation of /Cpeak and /Cnfw. With the above suitable modifications, our model works as long as the 
one-to-one correspondence of peaks and halos holds. We simply need to calibrate the mean scaling relation of /Cpeak.obs and 
ICpeak,h in absence of noise. 


APPENDIX C: FULL SKY RAY-TRACING SIMULATION 


Here we first summarize basic equations of the multiple-plane gravitational lensing algorithm, and then describe the ray¬ 
tracing method through the multiple-plane. For the former we largely follow [Pas fc Bode| ( |200^ , and for the latter we adopt 
one developed by Teyssier et al. (20091. 

Throughout this section, we work on the comoving coordinates. Thus p, x, and r(x) denote the comoving matter density, 
the radial comoving distance, and the comoving angular diameter distance, respectively. 


Cl Construction of the lensing potentials of multiple-plane 

The 3-dimensional light-cone matter distribution is composed of the multiple-layer of shells with a fixed width of 150 /i“^Mpc 
taken from the nested simulation boxes. The surface matter density field on a sphere for j-th shell is defined by 

Ai(0)=/ dx(p(0,x)-p)r■(x)^ (Cl) 

J shell 

where G denote the angular directions and p corresponds the mean matter density, and the integration is over the shell width. 
We set the lens-planes for each shell at the cone-volume weighted mean distance, Xj ~ 0.75(Xj,max “ Xymin)/(x|,max “ Xymin)) 
where Xi.max and Xi.min are the farthest and nearest radial distance to a shell, respectively. The convergence field for j-th 
shell is given by 


K\6) 
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(C2) 


where aj is the scale factor at the lens-plane Xj- We use the HEALPix ( |G6rski et al.|2005| ) scheme for pixelization of a sphere, 
and we make the best use of the HEALPix library. For each shell, we construct the projected mass density map from Y-body 
particles using Nearest Grid Point method (implemented with HEALPix subroutine vec2pix_ring). Using the volume (Kim) 
and total number of Y-body particles (Ypart) of the Y-body simulations, the total number of pixels (npix), the number of 
Y-body particles within i-th pixel (ripix,*) and its mean value (hpix), the convergence held is given by 


K^{e^) 
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(C3) 


Having the convergence held being ready, we expand it in spherical harmonics to have its coefhcients, K^, using the HEALPix 
subroutine map2alm. Then, the spherical harmonics coefhcients for the lensing potential can be obtained via 

for I / 0, (C4) 

and cpA ~ 0 for I = 0. This gives us the lensing potential held on a sphere, and its 1st and 2nd derivatives relate to the 
gravitational lensing dehection held (a^ = and the optical tidal matrix (U^^. = <()■’), respectively. Note that 

Vftj {i = 1,2) denotes the angular derivative. In an actual computation, we utilize the HEALPix subroutine alm2map_der. 
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C2 Light ray propagation 


Let us first describe the method to trace the ray trajectory using the multiple-plane algorithm, for which we basically follow 


one developed by Teyssier et al. (20091. A virtual observer is located at the center of the nested simulation boxes. Rays are 


traced backward from the observer point with the initial ray directions being set on HEALPix pixel centers. Thus the ray 
positions on the 1st (closest to the observer) is exactly at the HEALPix pixel centers. At each lens-plane, the ray directions are 
deflected according to aj, and the ray positions on the next lens-plane are computed using the method described in Appendix 
A of [Teyssier et ah] ( |2009| ). Note that ray positions on the j-th {j > 1) lens-plane are not exactly at the pixel center due to 
the lensing deflections, however the lensing fields (a^, are only computed at the pixel centers. In order to evaluate the 
lensing fields at an arbitrary position (0), we adopt the inverse distance weighted interpolation from nearest four pixel values; 


c4(0} = 


Efe=i WkaPdk -)■ e) 


e: 


Wfc 


(C5) 


where Wk = 1/|0 — 0fc|, and al{6k —>■ 0) is the deflection angle at the pixel center Ok but after parallel transporting to the 
ray position 0 (in order to take into account the change in the local (eg, e^) basis). In the actual computation, the parallel 
transport of the vector and tensor (for Ufj. described below) is implemented by the rotation of them by an angle between two 
coordinate bases at Ok and 0 (see Appendix C of |Becker|2013p . 

Having evaluated the ray positions at a lens-plane, we are able to compute the lensing magnification matrix (A^j.) using 
the recurrence relation ([Hilbert et al.|2009||Becker|2013||; 
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for j > 1, and note that in our notation, the lens-plane closest to the observer is j = 1. The optical tidal matrix, [//^, in 


the above relation is evaluated at the ray position on each lens-plane in the same interpolation scheme as eq. (C5|, and then 


again parallel transporting to the unperturbed ray position (i.e., the initial ray direction), because the observed magnification 
matrix should be evaluated in the local basis of the image position. 

We choose the source-plane at an arbitrary redshift Za (and thus the correspondence radial distance to the source-plane 
Xs), and evaluate the source position on the source-plane and the magnification matrix by using the above methods but 
replacing, e.g., Xi+i -t 


C3 Image positions of haloes 


The 3-dimensional light-cone distribution of dark matter haloes is generated in the same manner as for the matter distribution. 
The spatial position of a halo is converted into the angular position where the subscript “s” means the source position. 

We search for the corresponding image position in the following manner. First, we search for the nearest ray to the halo 
source position on the lens-plane of the shell where the halo is located. The displacement vector between the angular positions 
of the halo and the nearest ray is computed, A0s = —0^^. This vector is parallel transported to the image position of the 

nearest ray O’^p', and we denotes it by A0j Then the image position of the halo is given by A0/. The last step 

is valid if the difference in the lensing deflection angles between ray-trajectory to the halo and the nearest ray is very small. 
The statistical properties of differences in the lensing deflection angles between nearby two rays (the, so-called, the lensing 


excursion angle) were studies in Hamana & Mellier (20011; Hamana et al. (2005). They found that the root-mean-square 
(rms) value of the lensing excursion angles of rays for Za = 1 with the separation of 1 arcmin is ~ 1 arcsec. This value can be 
considered as the typical error in 0p^°. Considering the fact that the pixel scale of the current ray-tracing simulation is ~ 1 
arcmin, we may conclude that the above approximation is reasonably valid. However it should be noticed that for rays gone 
through a strong lensing region, the excursion angle can be much larger than the rms value, and thus may not be very 

accurate. There is room for improvement on this issue that we leave for future work. 


APPENDIX D: STATISTICAL PROPERTY OF PSEUDO-SPECTRUM ESTIMATORS 

In this appendix, we summarize the statistical property of pseudo-spectrum estimators. The pseudo-spectrum method is a 


powerful framework to construct the power spectrum of an underlying random held on limited sky (e.g., Hansen & Gorski 
2003 Efstathi^|2004 Brown, Castro &: Taylor|2005 ). 


Let us consider the two random helds in each direction in the sky: convergence held k{Q) and number density held of 
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lensing peaks p{^). These two fields would commonly be expanded in spherical harmonic as follows: 


(Dl) 


where yimiS^) represents the spherical harmonics and we can define ptm for the random field p(fl) similarly. The inverse 
transform is then given by 


— J' dn 


and the similar relation can be adopted for ptm- 

The effect of finite sky coverage for each field is characterized as 

k{Q) = W'"{Q)k{Q), 

p{n) = W’^{n)pin), 


(D2) 
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where and are the window function of sky masking for k and p, respectiveljQ Thus, the harmonic modes in presence 
of masked region is expressed as 


— ' m.' Xv T: 


where X = k, p and is defined by 

The estimators of power spectra on limited sky is defined by 

m 

where Xim = {kem,pim) and {XernXl^} represents the following set of power spectra: 




{k^rnP£m} 
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For the underlying field X = {k,p), we can define the power spectra P{i) as 
P{£) = lyXtraX\i^,)&u'^mm'■ 

Using Eqs. |D6[) and (|D9[), we can find the relation between P{£) and P{£) as follows: 
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The matrix Mu' represents the mode coupling effect due to masked region on power spectra, which is given by (in terms of 

P(^) = (P„„(^),Pp.(^),Ppp(€))^), 


Mu' — 


( 0 0 
0 0 
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(D13) 


0 Mf[, 


^ In practice, are not equal to This is because peaks of convergence field are defined by that of smoothed convergence map. 
When the area with mask is smoothed, there would exist ill-defined pixels due to the convolution between and a filter function 
for smoothing. Therefore, we need to remove the ill-defined pixels to find peaks. This procedure makes the effective sky coverage of 
smaller than that of W^. 
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where 
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and is given by 
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Hence, we can construct the estimator of P{£) from P{£) as 
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where P{i) is so-called psuedo-spectrum estimators. 

Next, we consider the covariance of the pseudo-spectrum estimators. The covariance of P{£) is defined by 
Coy[Pxy{£), PaInP')] = {PxY{t)PMN {(■')) — {PxYp)) {PmN {()) , 


(D19) 


where X, Y, M, N is set to be k or p. When the underlying held follows non-Gaussian and there exist no masked regions, the 
covariance can be expressed as 


Cov[Pxv(l^), PAf]v(^0]all-sky = 


21+1 


\PxM{t-)PYN{t') + PxN{t)PYM{t')\ + 


1 


2^ -f 1 2^' -b 1 


E (-^rmYrmM«'m'lV//^/)c,(D20) 


where the hrst term of the right-hand side in Eq. (D201 represents the Gaussian contribution to the covariance matrix and 


the second term corresponds to the contribution of four-point correlation function due to non-Gaussianity in the underlying 
held. On the other hand, in presence of masked region, the covariance of the pseudo-spectrum estimators is expressed as (see 
also, e.g., [Brown, Castro &: Taylor| ( [^05[ )) 

Cov[Pxv(^l), .P*rjv(^2)] = Cov[Pxv(^l), PAfiv(£2)]NG + 

e'+'^ 


where all the possible combinations of A, B, C and D are taken into account in Eq. (D211 and is given by 


y[XA,ND,MC,YB] _ _t_ 

“ (2£-b 1)(2P -b 1) 




(D22) 


The summation in Eq. (D22| is taken over all m, m', mi, m 2 . Here, denotes = 0, W'^'^ = W^, and 


y^rpp — y/p Eq. ( |D6| |. The non-Gaussian term in Eq. (D211 is dehned by 

1 

( 2 £i -b 


Cov[Pxv(^i),PAfiv(^ 2 )]NG = I l){2i2 + 1) 


MN\-1 
M' 




where the second summation in Eq. ( |D23 1 is over all mi, m', m'', £', I” {i = 1, 2) and the values of A, B, C and D. Eqs. (D211 
and (D231 clearly show that the complicated masked regions on sky would induce the additional mode-coupling of the 


covariance matrix of the pseudo-spectrum estimators. 
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